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This letter addresses basic questions concerning ferroelec- 
tric order in positionally disordered dipolar materials. Three 
models distinguished by dipole vectors which have one, two or 
three components are studied by computer simulation. Ran- 
domly frozen and dynamically disordered media are consid- 
ered. It is shown that ferroelectric order is possible in spatially 
random systems, but that its existence is very sensitive to the 
dipole vector dimensionality and the motion of the medium. 
A physical analysis of our results provides significant insight 
into the nature of ferroelectric transitions. 

PACS numbers: 64.70.Md, 77.80.-e, 82.20Wt 



Recently, spatially disordered dipolar materials have 
attracted considerable attention. Diluted lattices (jfl-fl, 
fluid phases ji^ and amorphous frozen ferrofluids 
have been examined experimentally pp jt^], theoretically 
P,[7[-pd]| and with computer simulations [ 3|-|5| . It has been 
shown that in the absence of long-range positional order, 
the strong spatial-oricntational coupling intrinsic to dipo- 
lar forces can lead to interesting new phase behavior. For 
example, frozen ferrofluids containing magnetic particles 
in a non-magnetic solvent exhibit magnetic irreversibil- 
ities reminiscent of randomly frustrated magnetic spin 
glasses [p^| . On the other hand, computer simulations of 
simple dipolar fluids clearly indicate the existence of a 
ferroelectric liquid crystal phase 

A simple interpretation of these observations might be 
as follows. In frozen ferrofluids, the quenched positional 
disorder creates random frustration and the system be- 
haves as a spin glass [jl2). The fluid systems |i|,|5| differ 
from the frozen case in that the strongly coupled transla- 
tional and rotational degrees of freedom are in full ther- 
mal equilibrium. This allows the development of short- 
range spatial correlations resembling those seen in the fer- 
roelectric tetragonal-I lattice and, consequently, fer- 
roelectric order develops in the liquid phase [Q. In view 
of these observations, and recalling that perfect crystals 
exhibit ferroelectric or antiferroelectric long-range order 
depending on the lattice symmetry [[U||lJ], one might 
argue that specific spatial correlations are required for 
ferroelectric order. 

In recent papers Zhang and Widom fill have put for- 
ward a mean field theory for spatially disordered dipolar 
systems. They argue that the long-range nature of the 
dipolar forces plays a key role in yielding ferroelectric 
order, and that this is not explicitly included in the sim- 



ple interpretation given above. More importantly, Zhang 
and Widom propose that, despite the strong frustration 
present in randomly frozen systems, long-range ferroelec- 
tric order is possible above a critical density. Their work 
implies that the spin-glass behavior observed in ferroflu- 
ids H results from the low concentration of magnetic par- 
ticles, whereas the ferroelectric liquid crystalline phase 
found in computer simulations of dipolar fluids Q|| arises 
because of the high particle density considered. In the 
present letter we examine the validity of this argument 
and address the general question "Can long-range fer- 
roelectric order spontaneously arise in a system without 
fined-tuned positional correlations?" 

We investigate the behavior of dense spatially disor- 
dered dipolar systems using constant temperature molec- 
ular dynamics (MD) and Monte Carlo (MC) simulations. 
Systems where the dipole vector has one, two and three 
components are considered. The first two of these are 
commonly referred to as the Ising and XY models and 
for notational simplicity we shall refer to the three com- 
ponent dipole as the XYZ model. In all cases the pair 
potential, u(12), is of the generic form 

u(12) = Ae{o/r) 12 +u DD {12) , 

where 4e(cr/r) 12 and 

u DD {12) = -3(>i • r)(>2 • r)/r 5 + ^ ■ fj, 2 /r 3 

are the soft-sphere and dipole-dipole interactions. The 
parameters e and a are the fundamental units of energy 
and length, \Xi is the dipole of particle i and r is the 
intermolecular vector. The long-range dipolar interac- 
tions were taken into account using periodic boundary 
conditions and the Ewald summation method assuming 
a perfectly conducting surrounding continuum ^],[^] . 
The existence of a ferroelectric phase can be detected by 
calculating the average polarization P per particle de- 
fined as P = l/N{\ X)<=i Aj ' where d is a unit vec- 
tor in the direction of the total instantaneous moment, 
M = YliLi Mii an d N is the number of particles in the 
system. 

It is convenient to characterize dipolar soft-sphere sys- 
tems by the reduced density, p* = Na 3 /V, the reduced 
temperature, T* = k^T/e, and the reduced dipole mo- 
ment, fi* — (fi 2 /eu 3 ) 1 / 2 , where V is the volume of the 
sample, T is the absolute temperature and fee is the 
Boltzmann constant. All results explicitly presented are 
for fi* — 4 and p* = 0.8. This density is well within the 
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range for which Zhang and Widom predict a ferroelec- 
tric phase Q. For (i* — 4 and p* = 0.8, Zhang and 
Widom predict a ferroelectric phase for the Ising case if 
T* < 35.2 and for the XYZ model if T* < 4.8. 

We first consider frozen systems. Suitable spatially 
disordered configurations were generated by carrying out 
an MD simulation of a soft-sphere fluid at T* = 10.5 
and p* — 0.8. Fluid-like configurations were then se- 
lected at random for dipolar rotational MD simulations. 
Following this approach we could obtain a frozen state 
at a much higher density than is possible from a ran- 
dom parking algorithm. Unfortunately, it is impossible 
to have a truly "random" and structureless spatial con- 
figuration (i.e., with the radial distribution function g(r) 
equal to one for r > a Jll[]) at this density. However, at 
T* = 10.5 the local structure in the soft-sphere fluid is 
weak and very short-ranged. 

Polarization results for randomly frozen systems are 
shown in Fig. 1 flf]. The XY and XYZ values were ob- 
tained with MD simulations. The discrete nature of the 
Ising model renders it inappropriate for MD so MC cal- 
culations were used. The average polarization obtained 
at the lowest temperature where equilibrium could be 
achieved, T^ in , is plotted vs 1/N. The values of T^ in 
are 10.0, 4.0, and 3.5 for the Ising, XY and XYZ models, 
respectively. Below these temperatures MD or MC runs 
for the same configuration started from perfectly ordered 
and disordered states (i.e., two replicas) did not converge 
to the same result within a reasonable computation time 
(i.e., about a week). Possibly with greater computational 
effort T* lin could be pushed a little lower, but the val- 
ues given above are well within the range where Zhang 
and Widom predict a ferroelectric phase. For the Ising 
model the equilibrated system at T* = 10.0 is nearly 
100% polarized. Moreover, for the Ising case the po- 
larization at T* lin exhibits very little number dependence 
and certainly does not appear to vanish in the thermody- 
namic limit. This, together with the plot of P vs T* and 
heat capacity calculations (see Fig. 3 below), strongly 
suggests that ferroelectric order develops spontaneously 
in the spatially disordered Ising system with the transi- 
tion occurring at T* w 25 for p* = 0.8. We note that 
this transition temperature is much lower than the value 
(T* = 35.2) predicted by Zhang and Widom. 

The situation for the XY and XYZ models is very dif- 
ferent. Although significant polarization was observed 
at finite temperatures, P decreases monotonously with 
increasing system size and appears to approach zero for 
an infinite system. Furthermore, the behavior of various 
spin-glass correlation functions and susceptibilities ]T^ ] 
suggests that both systems are entering an anisotropic 
spin-glass phase at nonzero temperature |20|, and that 
the observed polarization for the XY and XYZ models 
is due to a combination of short-range ferroelectric cor- 
relations and finite-size effects. Test calculations for the 
XYZ model using a denser frozen soft-sphere configura- 



tion (i.e., p* = 1.05 Q, T* = 10.5) also showed no 
long-range ferroelectric order. In brief, for the randomly 
frozen XY and XYZ models we find no evidence of a fer- 
roelectric state in the thermodynamic limit. This clearly 
disagrees with the calculations of Zhang and Widom 
which for the XYZ model predict a stable ferroelectric 
phase well within the temperature-density range consid- 
ered here. 

In order to gain further insight into the nature of fer- 
roelectric order (or the lack thereof) in spatially random 
systems, we consider a "toy model" where the transla- 
tional motion is completely decoupled from the dipolar 
interactions. The soft-sphere fluid acts as a "substrate" 
which moves at a fixed translational temperature inde- 
pendent of the embedded dipoles. The dipoles themselves 
interact and are equilibrated at a different rotational tem- 
perature. Of course, the "equilibrium" state achieved by 
the dipoles will depend on the underlying motion of the 
substrate. This model is similar in spirit to those used 
in recent studies of non-equilibrium phase transitions in 
magnetic systems subject to Levy flights [^l). It must be 
emphasized that this technique is presented only as a use- 
ful simulation tool and we do not imply any real physical 
mechanism for the decoupling. The moving substrate is a 
means of simulating dipolar systems in a dynamically ran- 
dom medium that lacks any specific spatial correlations 
which may favor ferroelectric ordering. The translational 
diffusion rate of the substrate can be controlled by the 
particle mass. Extrapolation to infinite mass should pro- 
vide information about the randomly frozen system. 

In Fig. 2, we have plotted P vs T* (rotational) for 
the XYZ model. Here, the decoupled substrate is a soft- 
sphere fluid again at p* = 0.8 and T* (translational) = 
10.5. It is convenient to introduce the reduced mass m* — 
m/m', where m! is a reference mass defined such that 
the reduced simulation timestep At* = (e/m'cr 2 ) 1 / 2 At = 
0.00125. Figure 2 includes results for m* = 1, 5 and 10. 
Spontaneous polarization develops for all systems and the 
temperature at which P begins to grow decreases with 
increasing mass. For example, from the P vs T* plot 
there appears to be a ferroelectric transition at T* ss 2 
for m* = 1. To verify that this is a real thermody- 
namic transition, we have calculated the Binder ratio 
111, SBindcr = (5/2) - (3/2)(|M| 4 )/(|M| 2 ) 2 , for systems 
with 64, 108 and 256 particles. A plot of ^Binder vs T* 
(see Fig. 2, inset) shows a clear crossing, and hence a 
transition, at T* — 1.9. Simulations were also carried 
out with m* = 20 but no significant orientational order 
was found above T* = 0.1. Very slow convergence pre- 
vented calculations at lower temperatures. 

We have also investigated dynamically disordered XY 
and Ising systems. The XY model behaves much like 
the XYZ system described above. For the Ising case, ro- 
tational dynamics cannot be used and a suitable Monte 
Carlo scheme which allowed the substrate to move in- 
dependent of the Ising dipoles was employed. P vs T* 
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results for m* = 1, m* = 5, and the randomly frozen 
system (m* = oo) are shown in Fig. 3. We see that for 
the Ising model the ordering behavior is essentially inde- 
pendent of mass and that the results for the dynamically 
disordered system with m* = 5 lie very close the those 
for the randomly frozen case. Heat capacities obtained 
by differentiating the average dipolar energy with respect 
to the rotational temperature are also shown in Fig. 3. 
The randomly frozen and m* = 5 results are very similar 
and indicate a phase transition at T* w 25. 

The dependence of the ferroelectric transition temper- 
ature on particle mass is shown in Fig. 4. The transi- 
tion temperatures were estimated from the heat capaci- 
ties and results are included for the XY and XYZ models. 
Results for the Ising system are not plotted because the 
transition temperature is essentially independent of the 
mass. We see that as the mass increases the transition 
temperature drops for both the XY and XYZ models. 
As noted earlier, for large masses and low temperatures 
convergence becomes prohibitively slow, but it seems rea- 
sonable to assume that the graph would simply continue 
with the transition temperature approaching zero in the 
infinite mass limit. This is consistent with the fact that 
we did not find a ferroelectric phase for randomly frozen 
systems at finite temperatures. 

The present and previous results can be interpreted as 
follows. It is useful to divide the local field, Ei oca i, ex- 
perienced by a particle in an infinite medium into two 
parts such that, Ei oca i = R + E, where R is a reac- 
tion field contribution ]j|[l5| and E is everything else. 
The reaction field arises from the long-range nature of 
the dipolar interactions; it is independent of the spatial 
correlations and favors ferroelectric order. The remain- 
ing contribution, E, is sensitive to positional correlations 
and may or may not favor ferroelectric order. Thus if 
R dominates a ferroelectric phase is to be expected, but 
if E is important the existence or nonexistence of a fer- 
roelectric phase will depend on the details of the spatial 
correlations. This simple picture allows us to rationalize 
the various systems considered. For fully coupled dipo- 
lar fluids [Q,H the short-range spatial correlations (and 
hence E) can adjust (i.e., become tetragonal- 1- like) in 
order to favor (or at least allow) a ferroelectric state. On 
the other hand, for randomly frozen systems the spatial 
correlations cannot adjust and at equilibrium E domi- 
nates R frustrating the formation of a ferroelectric phase 
except for the Ising case. Apparently, the discrete nature 
of the Ising model makes it much less susceptible to the 
development of disordering fields than are the XY and 
XYZ systems. The reaction field dominates in the Ising 
system and gives rise to a ferroelectric state. 

The behavior of the dynamically disordered systems 
can also be understood. If the translational diffusion of 
the substrate is sufficiently fast compared to the dipolar 
reorientational time, the extent of spatially dependent 
frustrating correlations is limited, R can prevail over E, 



and a ferroelectric phase is possible. As the mass of a 
substrate particle is decreased, the translational motion 
becomes faster and the above condition is met at higher 
and higher rotational temperatures. Thus the observed 
transition temperatures increase with decreasing mass. 
As the substrate particles become sufficiently light, the 
transition temperature is determined only by the reac- 
tion field and hence becomes independent of mass. In 
fact, m* = 1 gives essentially this limiting behavior and 
reducing the mass further has little effect on the transi- 
tion temperature. 

In conclusion, the answer to the question raised at the 
outset is a subtle one. Our results for the frozen Ising 
system and for the dynamically disordered XY and XYZ 
models clearly demonstrate that it is possible to have 
ferroelectric states without fine-tuned positional corre- 
lations. However, they also show that if a ferroelec- 
tric phase is to exist in a positionally random system, 
the long-range spatially-independent correlations arising 
through the reaction field must dominate the shorter- 
ranged position-sensitive correlations which generally act 
to frustrate ferroelectric order. For the Ising system 
the discrete nature of the interaction (i.e., there is just 
not much opportunity for favorable interactions among 
neighbors) limits the build up of short-range disordering 
correlations and there is a clear ferroelectric transition. 
For the dynamically random systems the orientational 
correlations opposing ferroelectric order are reduced by 
the motion of the substrate resulting in a ferroelectric 
phase. On the other hand, for the randomly frozen XY 
and XYZ models we find no indication of a ferroelec- 
tric phase at finite temperature. Rather, the disordering 
fields dominate and these systems appear to form nonfer- 
roelectric spin glasses at low temperature. Evidence for 
this is provided both by our direct simulations of frozen 
systems and by the extrapolation of the dynamically dis- 
ordered results to infinite mass. This conclusion must 
remain a little tentative because direct MD simulations 
at very low temperatures are not practical. Nevertheless, 
a ferroelectric phase in the randomly frozen XY and XYZ 
models seems highly unlikely. 
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FIG. 1. The polarization P at T^ in vs 100/iV for the 
randomly frozen Ising (circles), XY (squares) and XYZ (tri- 
angles). Results are included for 108 (XYZ only), 256, 392 
(XYZ only), 500 and 864 particles. 

FIG. 2. P vs T* (rotational) for dynamically random 
XYZ systems. The squares, triangles and circles are for 
m* = 1, 5 and 10, respectively. The error bars represent 
one estimated standard deviation. ^Binder vs T* (rotational) 
is shown in the inset for N = 64 (squares), 108 (triangles) 
and 256 (circles) particles. 

FIG. 3. P vs T* (rotational) for the Ising model. Re- 
sults are shown for dynamically random systems with m* = 1 
(squares) and m* = 5 (triangles) and for the randomly frozen 
case (circles). The heat capacities per particle, C v /N, are 
plotted vs T* (rotational) in the inset. 

FIG. 4. The mass dependence of the disordered to ferro- 
electric transition temperature Tp (rotational). The squares 
and triangles are results for the XY and XYZ models, re- 
spectively. The values of Tp were obtained from plots of the 
heat capacity, C v /N, vs T* (rotational) and a typical exam- 
ple (XYZ model, m* = 1) is shown in the inset. The error 
bars represent estimated uncertainties in the peak position. 
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